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ABSTRACT 



Context. One of the ways to determine the contribution of the dark halo to the gravitational potential of a galaxy is the study of non- 
circular (streaming) motions and the associated gas shocks in the bar region. These motions, determined by the potential in the inner 
parts, can break the disk-halo degeneracy. Here, two main fluid dynamical approaches have been chosen to model the non-circular 
motions in the bar region; a 2-D Eulerian grid code for an isothermal gas (FS2) and a 3-D smoothed particle hydrodynamic code 
(N-body/SPH) 

Aims. The aim of this paper is to compare and quantify the differences of the gas flows in rotating barred potential obtained using 
two different fluid dynamical approaches. We analyse the effect of using 2-D and a 3-D codes in the calculation of gas flow in barred 
galaxies and to which extend the results are affected by the code. To do this, we derive the velocity field and density maps for the 
mass model of NGC 4123 using a 3-D N-body/SPH code and compare the results to the previous 2-D Eulerian grid code results. 
Methods. Numerical modelling, 3-D N-body/SPH simulations 

Results. The global velocity field and the gas distribution is very similar in both models. The study shows that the position and 
strength of the shocks developed in the SPH simulations do not vary significantly compared to the results derived from the 2-D FS2 
code. The largest velocity difference across the shock is 20km s -1 between the 2-D and 3-D fluid dynamical models. 
Conclusions. The results obtained in the studies deriving the dark matter content of barred galaxies using the bar streaming motions 
and strength and position of shocks are robust to the fluid dynamical model used. The effect of 2-D and 3-D modelling can be neglected 
in this type of studies. 

Key words. Methods: numerical - Galaxies:kinematics and dynamics- Galaxies: structure - dark matter 



1. Introduction 

There have been several studies addressing the distribution of 
dark matter in galaxies using the non-axisymmetry of the poten- 
tial of ba rred galaxies : a) the fluid dynamics approach with the 
work bv lWeiner etail d2001al her eafte r WS WO 1) for the mod- 
elling of NGC 4123; iPerezI d2003l) and dPerez et al.Ll2004 here- 
after PFF04) for the study of the dark matter content of 5 barred 
galaxies (NGC 5505, NGC 7483, NG C 5728, and NGC 7267) ; 
and b) the sticky-particle approach by Rautiai nen et al. (2004) 



Salo et al 



1999) 



modelling the dynamics of ESO 566-24, and 
who modelled the Ha velocity field of IC 4214 to derive the halo 
contribution. 

The streaming motions and the associated gas shocks in the 
bar region are determined by the potential in the inner parts of 
galaxies. The velocity field of an axisymmetric galaxy does not 
allow to uniquely disentangle the contribution of the halo from 
that of the disk.The signatures of non-circular motions, however, 
can break the disk-halo degeneracy and be used to obtain the 
contribution of the dark halo to the potential. Common to all 
the modelling, the stellar potential is generated directly from the 
broad band galaxy images with some recipe to treat the verti- 
cal distribution and some assumptions about the mass-to-light 
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(M/L) of the stellar populations. The velocity fields obtained in 
this way are then compared to the observed kinematic informa- 
tion to determine if they can be reproduced by the models. 

The main result, common to all the modelling carried out up 
to now, is the fact that the gravitational field in the inner region 
is mostly provided by the stellar luminous component. The bar 
pattern speeds found by the different groups are consistent with 
fast rotators, and the best fitting M/L ratios obtained are compat- 
ible with M/L ratios derived from current population synthesis 
models. Although it is a very powerful method, only a handful 
of galaxies have been modelled in this way. 

Two main fluid dynamical approaches have been chosen for 
this type of modelling: an Eulerian grid code for an isother- 
mal gas (FS2) and a smoothed particle hydrodynamic code (N- 
body/SPH). SPH is a Lagrangian method for solving Euler's 
equation of motion. In this technique, the system is represented 
by a set of particles and the gas properties are calculated by a 
weighted average over the neighbouring particles. On the other 
hand, the Eulerian grid is fixed in space and time, the grid nodes 
and cells remain spatially fixed while the materials flow through 
the mesh. In the modelling carried out by WSW01, the scale 
height is introduced as a smoothing length for the potential while 
the approach chosen by PFF04 uses a full 3-D calculation of the 
forces. For future work and as a consistency check, it is impor- 
tant to understand whether the results obtained are dependent a) 
on the dynamical code used and b) 2D versus 3-D. No significant 
differences are expected between the SPH and the Eulerian grid 
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approach. Previous simulati ons (Englmaier & Gerhard fl997t 
iPatsis & Athanassoulail2000l hereaft er PAOO), running 2-D SPH 
under the same conditions as used by I Athanassoulal 11992) with 
the FS2 Eulerian grid code indicated that the main features, such 
as the strength of the shocks in the bar are quite similar and small 
differences arise due to the statistical nature of SPH or to numer- 
ical and artificial viscosities used in the different codes. Both 
studies pointed out the importance of certain parameters used in 
the simulations for the outcome of the modelling, such as the 
sound speed and the way the non-axisymmetric part of the po- 
tential is introduced. The use of a 2-D versus a 3-D code might 
introduce large variations in the radial forces due to the differ- 
ent treatment of the vertical forces, which in turn might vary 
the outcome of the analysis and the conclusions derived from 
this methodology. In this paper, we investigate the effect of 2-D 
vs 3-D numerical modelling in the gas flows of rigidly rotating 
bar potentials. In order to test this effect, we have modelled the 
gas flow in NGC 4123 using the mass model and images pro- 
vided by B.J. Weiner but using the 3-D N-body/SPH code used 
in PFF04. The parameters chosen for the SPH simulations have 
been taken from WSW01. In future work we will address the 
comparison between the sticky particle codes and the fluid dy- 
namical methods. The modelling is presented in Section [2] the 
results and comparison are presented in Section [3] Finally, dis- 
cussion and conclusions are presented in Section|4] 

2. Modelling 

The methodology for this study is simple, we derive the mass 
model for NGC 4123 in a way similar to that of WSW01. We 
then run the SPH models to obtain the gas distribution and ve- 
locity field, using the best fit parameters (M/L and bar pattern 
speed) in WSW01. And finally, the gas distribution and the ve- 
locity field are compared to the results presented in WSW01. In 
this section we will present a brief summary of both codes and 
initial conditions, emphasising the SPH carried out in this work 
and explaining the way in which we derive the mass distribution. 

The FS2 code is a second-order, flux splitting, Eulerian grid 
code for an isothermal gas in an imposed gra vitational potential, 
originally written by G.D. Ivan Albadal d 1985b . The sound speed 
of the gas is 8km s _1 .The code does not include self-gravity of 
the gas. The barred potential rotates at a fixed pattern speed Q#, 
and the grid rotates with the bar. In WSW01, 0.1 Gyr are al- 
lowed to fully grow the bar in order to steadily adjust the flow. 
The initial gas surface density is set to be 1OM pc~ 2 inside 
a radius of 8 kpc; outside that radius it falls off exponentially. 
A sech 2 (z/zo), where zo is the scale height, is assumed for the 
vertical distribution. With this distribution, the accelerations are 
calculated in the midplane using a Fourier transform method to 
conv olve a Green's function with the surface brightness distribu- 
tion |Hockne3([l965|). The scale height is effectively a smoothing 
length for the potential. For a complete explanation of the code 
and the initial conditions refer to WSW01. 

The 3-D N-body and hydro code was initially developed 
by the Geneva Ob servatory galactic dynamics group for spi - 
ral galaxy studies dPfenniger & Friedlil Il993l; IFuxL 1 1 997L 1 1 999h . 
The, initially self-consistent, code was modified to use a fixed ro- 
tating potential. The stellar potential is fixed using the observed 
light distribution. The potential is calculated using a particle- 
mesh technique with cylindrical-polar grid and the short range 
forces are softened using a variable homogeneous ellipsoid ker- 
nel with principal axes matched to the local grid resolution. The 
pressure forces and viscous fo rces are der i ved by 3- D Smooth 
Particle Hydrodynamics, SPH dBenzj|1990tlFuxlll999l) . The gas 



is taken to be isothermal with a sound speed of 10km s _1 and an 
adiabatic index of 5/3, corresponding t o the ato mic hydrogen. 
For mo re de t ails ab out the code, refer to lFuxl ( 1 19991) or his Ph.D. 
Thesis. IFuxI (119971) . For details about the code in the form used 
for this work refer to PFF04. 

The grid size used is 95 cells in the radial, 96 in the azimuthal 
and 1214 (including doubling up) in the vertical directions. The 
vertical resolution is set to 0.05 times the scale-height adopted 
for the luminous mass distribution. The number of gas particles 
used is 300 000. The barred potential rotates at fixed pattern 
speed Q.b and there is no self-gravity. These simulations were 
run; first, using the Alpha Server Sc system at the Australian 
National University super computer facility; and later, the Cray 
SVle super computer at Groningen University. 

For the initial conditions see PFF04. The initial gas density 
in the simulations consists of one component. The radial distri- 
bution for the gas is a Beta function with a standard deviation set 
to the scale-length of the visible disk and radially vanishing at a 
distance 4 times this scale-length. The scale-length assumed for 
NGC 4123 is 3.2 kpc, as derived from the /-band luminosity dis- 
tribution in lWeiner et al.l(l2001bl) . The vertical distribution of the 
gas is generated directly by solving the hydrostatic equilibrium 
equation for an isothermal gas. The gas particles have pure circu- 
lar motion with cylindrical rotation and zero velocity dispersion, 
since the effective dispersion is taken into account in the pressure 
component of the SPH. The bar is chosen to grow linearly during 
three bar rotations, so that the gas flow can steadily adjust to the 
forcing bar without requiring much CPU time. In previous sim- 
ulations (PFF04) different onset times were checked to ensure 
that at the chosen onset time the particle flow adjusted steadily, 
finally setting the onset time to 3 bar rotations. 

We have used the /-band light distribution, and exponen- 
tial vertical profile with a scale-height of 200 pc to derive the 
potential in a similar fashion as in WSW01. They constructed 
a bi-symmetric image by rotating the image 180° and averag- 
ing the rotated and original images. The central point source 
was removed and modelled as a spherical Gaussian distribution 
p(r) oc exp(-r 2 /2r 2 ) with a scale radius (r c ) of 200 pc, refer 
to WSW01 for details about the mass distribution. Since the 
model used for comparison is the maximum disk model, no dark 
matter halo profile was included in the derivation of the mass. 
Because the potential in the bar region is dominated by the stel- 
lar component we do not expect a significant difference in the 
results when excluding this component. The model parameters 
used were M/L of 2.25 and pattern speed of 20 km s^'kpc -1 , 
corresponding to a /? cor /Rbar = 1-35, where R co , is the corotation 
radius and Rb ar is the bar semi-major axis. These parameters cor- 
respond to the best-fit model in WSW01. To derive the velocity 
field, we used the parameters derived by IWeiner et alJ (2001b) 
from their kinematic data; i.e. i = 45 ° and major axis of the 
galaxy at 57 ° North through West. 

Both codes were run without gas self-gravity. Low resolution 
simulations have sh own that gas self-gravity i s irrel evant for the 
results in this work dPerezl 120031: iPerez et al.L |2004|) . 

The gas density distribution, for the SPH modelling, clearly 
settles into a stationary configuration after 3 bar rotations, a nu- 
clear ring develops perpendicular to the bar, associated to the X2 
family of orbits . From the frequency plot of the axisymmetric 
case, we see that NGC4123 possesses two ILR for the derived 
pa ttern speed. The pres ence of the ILR was already suggested 
in IWeiner et al.1 (12001 bl) from the offset of the dust lanes in the 
colour images of NGC 4123. 
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Fig. 1. On the left panel; a) the velocity field contours in the bar region, projected into the plane of the sky, using i = 45 and 
major axis of the galaxy at 57 ° northwest from WSW01. Superposed are the isovelocity contours of the projected velocity field as 
derived in this work, with levels corresponding to intervals of 30 km s _1 . On the right panel; b) density distribution from WSW01, 
superposed are the isodensity contours from the particle distribution of the SPH for the mass distribution of NGC 4123 derived from 
the /-band after 4 bar rotations. The non-axisymmteric component is fully grown after three bar rotations and the bar patter speed is 
rotating clockwise in the inertial frame. Notice the agreement in the regions of high density within the bar region. 



3. Comparison and results 

After three bar rotations, when the non-axisymmetric component 
is fully grown, the gas distribution settles into a steady configu- 
ration until the end of the simulation, six bar rotations. The gas 
distribution at the end of the simulation is shown in Fig. 1. The 
region chosen for a comparison between the two models is the 
region used by WSW01 in their study (see Fig.[T]for the location 
of the comparison regions). They chose the bar region because 
is the area where the non-axisymmetric motions are strong and 
excluded the central region because the velocity gradients are 
very high and probably unresolved. Their dynamical resolution 
is worst at the center because there are few cells across the scale 
of a streamline. Therefore, for a meaningful comparison between 
the models, we will restrict the comparison to the regions dis- 
cussed in WS WO 1. 

Previous work (PA00) in a study of the effect of hydrody- 
namical and numerical parameters in the simulations of gas flow 
in barred galaxies, noticed that differences in the sound speed 
would change the gas distribution morphology and the inflow 
rate of the gas. The two models under comparison here have sim- 
ilar sound speeds (10km s _1 for the SPH models, and 8km s~ ! for 
the grid based models) and therefore, we do not expect differ- 
ences arising from this effect. 

Other parameter that could cause a difference between the 
two models is the way the non-axisymmetric component is in- 
troduced, PA00 noticed that introducing the bar abruptly results 
in a different morphology when compared to the case where the 
bar is slowly introduced. In both models analised in here, the bar 



is introduced slowly, allowing the flow to steadily adjust to the 
growth of the bar. 

The spatial gas resolution in the two codes is very different 
and difficult to compare. The gas resolution in the SPH is deter- 
mined by the particle smoothing length, in this case, to ensure 
good resolution in high density regions, the smoothing length is 
assigned to each particle such that their number of neighbouring 
particles always remains close to a fixed number. The average 
gas smoothing length is close to the vertical gravitational reso- 
lution, which for NGC 4123 is 20 pc. However, in the central 
parts, where the particle density is high, the resolution reaches 
a few parsecs. The resolution of the grid based code is 84.7 
pc 2 , using a grid having 256x512 cells. It is interesting to no- 
tice that while the gas distribution in the SPH code develops a 
central ring, associated to the ILR, in the results obtained with 
the Eulerian code the gas density distribution does not present 
this central feature; although, the underlying potential and the 
sound speed are similar in both models. Because the mass dis- 
tribution and the bar pattern speed are similar, the difference in 
the presence or not of the central ring is possible due to differ- 
ences between the two codes. PA00 showed the gas distribution 
for a series of 2-D SPH models with different gas sound speeds 
and compared them to simulations carried out with a 2-D sec- 
ond order flux-splitting code similar to the one used by WSW01 . 
The gas response for a sound speed of lOkms -1 was similar in 
both cases, both cases developing a nuclear ring. However, for a 
sound speed of 15km s _1 , the gas density response is different in 
both codes. The gas density in the SPH code develops a central 
X2 ring, while the gas distribution in the eulerian code gives no 
central ring. This should be further investigated using 2-D (or 



4 



I. Perez: Comparison between 2-D and 3-D codes in dynamical simulations of gas flow in barred galaxies 






Fig. 2. This figure shows two cuts across the high velocity gradient region for the deprojected velocity maps of both models. The 
top panels show the projected velocity. The position corresponds to the velocity going from the leading side towards the trailing 
side of the shock region in a 30 pixel range {approx 23.4 arcsec), with a cut a a distance along the shock of 30 pixels from the 
centre (left panel) and a distance of 35 pixels from the centre (right panel). The middle panels correspond to the residuals of both 
line-of-sight velocity distributions. The bottom panels shows the position of the cuts (two parallel black lines) shown in the above 
panels perpendicular to the shock region on top of the general gas distribution. 
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3-D) consistently. The comparison we carry out throughout this 
paper is done outside this central region. Future high resolution 
observations of the gas distribution in the central parts, would 
show if the velocity fields and the density distribution generated 
in our modelling is a good representation of the central region 
in NGC 4123 or whether we are making wrong assumptions on 
the vertical thickness of the bar or the dynamical approach used 
here is not appropriate for the problem of modelling the central 
regions. 

The general characteristics of the velocity fields are the fol- 
lowing; the offset shocks in the bar are roughly east-west and 
nearly perpendicular to the bar there is a large velocity gradient. 
In Fig.QJ) one can see that, within the comparison region, i.e. the 
regions associated to the dust lanes the agreement between both 
models is very good. These regions are associated to shocks and 
to regions of high velocity gradients. 

The velocity field and the corresponding velocity gradients 
are similar in both models as shown in Fig.QJ. The shape of the 
density maxima (accompanied by the hi gh velocity gr a dients ) 
puts strong constrains to the potential dAthanassoulal Q992). 
The angle between the high velocity gradient region and the 
bar is the same in both simulations and the shape of the shock 
regions in both models is similar, supporting the results from 
both models. There seems to be a slight displacement of the 
peak position of the large gradient region between both mod- 
els (« 1.5 arcsec), in the SPH the regions of highest density are 
broader than in the FS2 code possibly due to the fact that the 
shock is slightly less well resolved in the SPH model. This dif- 
feren ce between SPH and the grid codes has been already no- 
ticed Engl maier & Gerhard ( 1997). They obtained broader den- 
sity peaks in the shock regions in 2-D SPH modelling compared 
to 2-D grid models. 

To compare more quantitatively the two models we have 
made cross section profiles of this region of high velocity gra- 
dients (see Fig. |2). To facilitate the comparison of the veloc- 
ity jumps in the shock region, the velocity distribution has been 
shifted by two pixels (1.56 arcsec) because this is the offset 
found between both models in the peak position of the gradi- 
ent. The cut positions have been chosen to be within the com- 
parison region. One has to keep in mind the statistical nature 
of SPH when comparing the cuts. This causes some particles to 
have very different values from the general distribution. The ve- 
locity jumps are very similar for both models with a maximum 
velocity difference between the two models across this region of 
20 kms 1 . The amplitude of the velocity jump depends on the 
mass distribution in the bar region, wi th thinner and more mas- 
sive bars having larger velocity jumps (lAthanassou kl ll992l) . and 
of course, in this case, also on the position of the cut across the 
shock. This is what gives rise to the difference between the right 
and the left panel in Fig. [2] Overall, the shape and the amplitude 
of the gas velocity cross section profiles between the two models 
agree very well. 

A change in the parameters such as the pattern speed and 
the M/L would have caused the velocity field and the shock 
strengthjocat i on an d shape to be signifi c antly different (se e 
lAthanassoulal ( 1 1992b : IWeiner et all d200Tah ; iPerez et all d2004h ) 
to the results obtained here. 

There is a discrepancy outside the bar region, and the posi- 
tion and shape of the spiral arms do not agree between the two 
models. In the Eulerian code the arms are narrower than in the 
SPH code; they seem to circumscribe the arms formed in the 
SPH code. In any case, the region of the spiral arms is discarded 
in the comparison with the observations, both groups assume a 



single pattern speed which is probable not the case for real galax- 
ies. 

4. Conclusions 

The position and strength of the shocks in barred galaxies de- 
rived from dynamical models is robust to the chosen model. 
Here, we have compared 2-D to a 3-D code and despite the dif- 
ferences in the two codes and the small differences in the mass 
model, the resulting velocity and density maps for both models 
is very similar. However, for similar hydrodynamical parame- 
ters, the SPH code develops a central ring, associated to the X2 
family of orbits while the gas distribution in the FS2 model does 
not present this central ring. This difference seems to be due to 
differences between the two codes and not to the difference on 
the treatment of the vertical forces. In any case, the position and 
the strength of the shocks in both models remains very similar. 
We can therefore conclude that the methodology used to derive 
the halo contribution to the potential using dynamical modelling 
of the gas flow in the bar region of spiral galaxies is robust to the 
codes used and to wether they are 2 or 3-D codes. 
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